#include "euler.hh"
#include <math.h>
namespace Dune {
  namespace Initial {
    Euler::Euler(Dune::ParameterTree& ptree){
      set_initcond_para("dim",1.0);
      set_initcond_para("sod_shocktube",0.0);
      set_initcond_para("blastwave",0.0);
      set_initcond_para("two_shockwave_1d",0.0);
      // 
      set_initcond_para("cen1",0.0);
      set_initcond_para("cen2",0.0);
      set_initcond_para("cen3",0.0);
      set_initcond_para("slope1",0.0);
      set_initcond_para("slope2",0.0);
      set_initcond_para("slope3",0.0);  
      set_initcond_para("radius",0.0);
      //
      set_initcond_para("rhoI",0.0);
      set_initcond_para("rhoII",0.0);    
      set_initcond_para("mom1I",0.0);
      set_initcond_para("mom1II",0.0);
      //
      set_initcond_para("eI",0.0);
      set_initcond_para("eII",0.0);
      set_initcond_para("gamma",5.0/3.0);
      //
      // std::cout << "Initial Euler initialization" << std::endl;
      set_initcond_paras(ptree); 
    }

    Euler::~Euler(){
    }

    void Euler::set_initcond_paras(Dune::ParameterTree& ptree){
      std::cout << "set_initcond_paras" << std::endl;
      para["sod_shocktube"]=ptree.get("initial.sod_shocktube",0.0);
      para["blastwave"]=ptree.get("initial.blastwave",0.0);
      para["sod_two_shockwave_1d"]=ptree.get("initial.two_shockwave_1d",0.0);
      // 
      para["dim"]=ptree.get("grid.dimension",0.0);
      para["cen1"]=ptree.get("initial.cen1",0.0);
      para["cen2"]=ptree.get("initial.cen2",0.0);
      para["cen3"]=ptree.get("initial.cen3",0.0);
      para["slope1"]=ptree.get("initial.slope1",0.0);
      para["slope2"]=ptree.get("initial.slope2",0.0);
      para["slope3"]=ptree.get("initial.slope3",0.0);
      para["radius"]=ptree.get("initial.radius",0.0);
      para["rhoI"]=ptree.get("initial.rhoI",0.0);
      para["rhoII"]=ptree.get("initial.rhoII",0.0);
      para["mom1I"]=ptree.get("initial.mom1I",0.0);
      para["mom1II"]=ptree.get("initial.mom1II",0.0);
      para["eI"]=ptree.get("initial.eI",0.0);
      para["eII"]=ptree.get("initial.eII",0.0);
      para["gamma"]=ptree.get("initial.gamma",0.0);
    }

    void Euler::initial_condition(double* location, double* holder){
      // std::cout << "choice initial condtion" << std::endl;
      if(para["sod_shocktube"]==1.0)
	sod_shocktube(location, holder);
      if(para["blastwave"]==1.0)
	blastwave(location, holder);
      if(para["two_shockwave_1d"]==1.0)
        two_shockwave_1d(location, holder);
    }

    // 1-dimension: Figure 14.3 of book: Leveque R. J.
    void Euler::sod_shocktube(double* location, double* holder){
      std::cout << "initial condtion: perform sod_shocktube" << std::endl;
      int dim = (int)para["dim"];
      double cen[dim], slope[dim];  
      double dotwithnormal=0;
      double tmp=0.0;
      para["cen1"]=0.5;
      para["cen2"]=0.5;
      para["cen3"]=0.5;
      para["slope1"]=1.0;
      para["slope2"]=1.0;
      para["slope3"]=1.0;    
      para["rhoI"]=3.0;
      para["rhoII"]=1.0;    
      para["pI"]=3.0;
      para["pII"]=1.0;
      if(dim>=1){
        cen[0]=para["cen1"];
        slope[0]=para["slope1"];
      }
      if(dim>=2){
        cen[1]=para["cen2"];
        slope[1]=para["slope2"];
      }
      if(dim==3){
        cen[2]=para["cen3"];
        slope[2]=para["slope3"];    
      }
      for(int i=0; i<dim;i++)
        dotwithnormal=dotwithnormal+slope[i]*(location[i]-cen[i]);
      holder[0] = dotwithnormal>0? para["rhoI"] : para["rhoII"];
      for(int i=1; i<=dim; i++)
        holder[i] = 0.0;
      holder[dim+1] = dotwithnormal>0? para["pI"]/(para["gamma"]-1.0) : 
        para["pII"]/(para["gamma"]-1.0);  
      return;
    }

    // Figure 14.4 of book: Leveque R. J.
    void Euler::blastwave(double* location, double* holder){
      int dim = (int)para["dim"];
      para["rhoI"]=1.0;
      para["pI"]=1.5;
      para["pII"]=1.0;
      para["radius"]=0.25;
      double distance=0;
//      std::cout << "use blastwave initializer" << std::endl;
//      std::cout << "rhoI: " << para["rhoI"] << std::endl;
//      std::cout << "pI: " << para["pI"] << std::endl;
//      std::cout << "pII: " << para["pII"] << std::endl;
//      std::cout << "radius: " << para["radius"] << std::endl;
      //
      holder[0] = para["rhoI"];
      for(int i=1; i<=dim; i++)
        holder[i] = 0.0;
      for(int i=0; i<dim; i++){
        distance+=std::pow(location[i]-0.5,2);
      }
      distance=std::sqrt(distance);
      holder[dim+1] = distance<=para["radius"] ? 
        para["pI"]/(para["gamma"]-1.0) : para["pII"]/(para["gamma"]-1.0);
      return;
    }

    void Euler::two_shockwave_1d(double* location, double* holder){
      para["cen1"]=0.3;  
      para["rhoI"]=1.0;
      para["rhoII"]=2.0;    
      para["mom1I"]=1.0*3.0;
      para["mom1II"]=2.0*1.0; 
      para["pI"]=1.0;
      para["pII"]=1.0;
      holder[0] = (location[0]<para["cen1"]) ? para["rhoI"] : para["rhoII"];
      holder[1] = (location[0]<para["cen1"]) ? para["mom1I"] : para["mom1II"];
      holder[2] = (location[0]<para["cen1"]) ? 
        std::pow(para["momI"],2)/para["rhoI"]/2.0+para["pI"]/(para["gamma"]-1.0) : 
        std::pow(para["momII"],2)/para["rhoII"]/2.0+para["pII"]/(para["gamma"]-1.0);
      return;
   }
 } // Initial
}  // Dune
